#library(ncdf4)
library(ggplot2)
require(data.table)
library(dplyr)
library(plyr)

#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\...\\NFood_replication\\"


#=========Supplement Figure S3 for sensitivity analysis with interaction terms

load(paste0(wd1,"/data/RData/Fig.ED4.MLR-pergrid-rad-ppt.log.Rdata")) #coefficients when log change of RADD,RADI,P,RH are predictors
coeff.i.sai <- coeff.i.samples3[Scenario=="SAI - RCP85"]
coeff.r.sai <- coeff.r.samples33[Scenario=="SAI - RCP85"] 
coeff.r.sai <- coeff.r.samples333[Scenario=="SAI - RCP85"]

#interannual variability under RCP85 (absolute change): -1sd T ~ -1.5 [°C], +1sd RH ~ +3.6 unit RH [%], -1sd P ~ -22 P [mm/month], -1sd RD ~ -5.5 [W/m^2], +1sd RI ~ 1.4 [W/m^2]
#interannual variability of SAI-RCP85: -1sd T ~ -1.06 [°C], +1sd RH ~ 4.4 [%] or 0.08 [log], -1sd P ~ -0.31 P [log], -1sd RD ~ -0.096 [log], +1sd RI ~ 0.08 [log]
dif.8crop.sai <- copy(dif.8crop.s[Scenario=="SAI - RCP85"])
dif.8crop.sai <- dif.8crop.sai[, tsa.sd:=sd(tsa.dif), by=.(Crop,Irrig,Sample)]  #absolute terms for T and RH
dif.8crop.sai <- dif.8crop.sai[, rh.sd:=sd(rh.dif), by=.(Crop,Irrig,Sample)]
dif.8crop.sai <- dif.8crop.sai[, radd.sd:=sd(radd.log), by=.(Crop,Irrig,Sample)] #log terms for RD, RI, P
dif.8crop.sai <- dif.8crop.sai[, radi.sd:=sd(radi.log), by=.(Crop,Irrig,Sample)] #because they have log linear relation with yield.dif (log)
dif.8crop.sai <- dif.8crop.sai[, ppt.sd:=sd(ppt.log), by=.(Crop,Irrig,Sample)]
dif.8crop.sai[,.(tsa.sd=weighted.mean(tsa.sd, area.wt, na.rm = T), rh.sd=weighted.mean(rh.sd, area.wt, na.rm = T), 
               ppt.sd=weighted.mean(ppt.sd, area.wt, na.rm = T), radd.sd=weighted.mean(radd.sd, area.wt, na.rm = T),
               radi.sd=weighted.mean(radi.sd, area.wt, na.rm = T))]
#use absolute changes in T and RH to estimate the effects of "-1 [°C] T", "+5 [%] RH", but log changes in RD, RI and P (they have log linear relationship with yield.dif)
#also show +- sd in all variables to estimate the effects of "-1sd T", "-1sd RD", "+1sd RI", "-1sd P", "+1sd RH"

#unit changes should be named c("-1 [°C] T", "-10% RD", "+10% RI", "-10% P", "+5 [%] RH"), and c("-1sd T", "-1sd RD", "+1sd RI", "-1sd P", "+1sd RH") respectively
#either apply the same unit changes to all grid crop grid cells or use crop-grid-specific sd 
dif.unit <- dif.8crop.sai[,.(tsa.unit1=-1, rh.unit1=5, radd.unit1=log(1-0.1), radi.unit1=log(1+0.1), ppt.unit1=log(1-0.1), 
                          tsa.unit2=-tsa.sd, ppt.unit2=-ppt.sd, rh.unit2=rh.sd, radd.unit2=-radd.sd, radi.unit2=radi.sd, #alternative: 1sd of interannual variability of SAI-RCP85 as unit changes 
                          #tsa.unit2=-weighted.mean(tsa.sd, area.wt, na.rm = T), radd.unit2=-weighted.mean(radd.sd, area.wt, na.rm = T),radi.unit2=weighted.mean(radi.sd, area.wt, na.rm = T),
                          #rh.unit2=weighted.mean(rh.sd, area.wt, na.rm = T), ppt.unit2=-weighted.mean(ppt.sd, area.wt, na.rm = T),
                          yield0, area.wt,Year,Crop,Irrig,Sample)]

sens.i.sai <- coeff.i.sai[dif.unit[Irrig==TRUE],                                    
                          .(Intercept, T.log1=T*i.tsa.unit1, RD.log1=RD*i.radd.unit1, RI.log1=RI*i.radi.unit1, RH.log1=0, P.log1=0, 
                            T.log2=T*i.tsa.unit2, RD.log2=RD*i.radd.unit2, RI.log2=RI*i.radi.unit2, RH.log2=0, P.log2=0, 
                            yield0=i.yield0, area=i.area.wt, Year=i.Year), by=.EACHI, on=.(Crop,Sample)]
sens.r.sai <- coeff.r.sai[dif.unit[Irrig==FALSE], 
                          .(Intercept, T.log1=T*i.tsa.unit1, RD.log1=RD*i.radd.unit1, RI.log1=RI*i.radi.unit1, P.log1=P*i.ppt.unit1, RH.log1=RH*i.rh.unit1, 
                            T.log2=T*i.tsa.unit2, RD.log2=RD*i.radd.unit2, RI.log2=RI*i.radi.unit2, P.log2=P*i.ppt.unit2, RH.log2=RH*i.rh.unit2, 
                            TRH.log1=TRH*i.rh.unit1*i.tsa.unit1, TRH.log2=TRH*i.rh.unit2*i.tsa.unit2, 
                            #RDRI.log1=RDRI*i.radd.unit1*i.radi.unit1, RDRI.log2=RDRI*i.radd.unit2*i.radi.unit2,
                            RHP.log1=RHP*i.rh.unit1*i.ppt.unit1, RHP.log2=RHP*i.rh.unit2*i.ppt.unit2,
                            yield0=i.yield0,area=i.area.wt, Year=i.Year), by=.EACHI, on=.(Crop,Sample)]
sens.i.sai  <- sens.i.sai[,Tot.log1:=Intercept+T.log1+RD.log1+RI.log1] 
sens.r.sai  <- sens.r.sai[,Tot.log1:=Intercept+T.log1+RD.log1+RI.log1+P.log1+RH.log1+TRH.log1+RHP.log1]
sens.i.sai  <- sens.i.sai[,Tot.log2:=Intercept+T.log2+RD.log2+RI.log2] 
sens.r.sai  <- sens.r.sai[,Tot.log2:=Intercept+T.log2+RD.log2+RI.log2+P.log2+RH.log2+TRH.log2+RHP.log2]

sens.i.sai$Irrig <- TRUE
sens.r.sai$Irrig <- FALSE
sens.sai <- rbind(sens.i.sai, sens.r.sai, fill=TRUE)

sens.sai.wm <- sens.sai[, .(T.log1=weighted.mean(T.log1, area, na.rm = T), RD.log1=weighted.mean(RD.log1, area, na.rm = T),RI.log1=weighted.mean(RI.log1, area, na.rm = T),
                           RH.log1=weighted.mean(RH.log1, area, na.rm = T), P.log1=weighted.mean(P.log1, area, na.rm = T), Tot.log1=weighted.mean(Tot.log1, area, na.rm = T), 
                           T.log2=weighted.mean(T.log2, area, na.rm = T), RD.log2=weighted.mean(RD.log2, area, na.rm = T), RI.log2=weighted.mean(RI.log2, area, na.rm = T),
                           RH.log2=weighted.mean(RH.log2, area, na.rm = T), P.log2=weighted.mean(P.log2, area, na.rm = T), Tot.log2=weighted.mean(Tot.log2, area, na.rm = T),
                           R.log1=weighted.mean(RD.log1+RI.log1, area, na.rm = T), R.log2=weighted.mean(RD.log2+RI.log2, area, na.rm = T), 
                           TRH.log1=weighted.mean(TRH.log1, area, na.rm = T), TRH.log2=weighted.mean(TRH.log2, area, na.rm = T), 
                           #RDRI.log1=weighted.mean(RDRI.log1, area, na.rm = T), RDRI.log2=weighted.mean(RDRI.log2, area, na.rm = T), 
                           RHP.log1=weighted.mean(RHP.log1, area, na.rm = T), RHP.log2=weighted.mean(RHP.log2, area, na.rm = T), 
                           yield0=weighted.mean(yield0, area, na.rm = T), area=sum(area, na.rm=T)),
                       by =.(Crop,Irrig,Year)]
sens.sai.wm$Level <- "Mean"

rm(dif.8crop.sai, dif.unit, sens.i.sai, sens.r.sai)


#####################load resampled data for confidence intervals######################
#see "Data.ED.Fig.4.R": do resampling to get confidence interval of sensitivity effect of each variable

load(paste0(wd1,"/data/RData/Fig.ED4.MLR-pergrid-sensitivity-SAI-Bootstrap.Rdata")) #include T*RH and RH*P
str(sens.log.boots)
unique(sens.log.boots$Year)

#25-year mean effect of 2075-2095 (same sensitivity across years but crop distribution changes, thus global average changes)
sens.sai.25yr <- sens.log.boots[Year>=2075 & Year<=2099, lapply(.SD, mean, na.rm=TRUE), by =.(Crop,Irrig,Level)] 

#translate log to % effect
sens.sai.pct <- sens.sai.25yr[, .( T.eff1=(exp(T.log1)-1)*100, RD.eff1=(exp(RD.log1)-1)*100, RI.eff1=(exp(RI.log1)-1)*100, P.eff1=(exp(P.log1)-1)*100, RH.eff1=(exp(RH.log1)-1)*100,
                                   T.eff2=(exp(T.log2)-1)*100, RD.eff2=(exp(RD.log2)-1)*100, RI.eff2=(exp(RI.log2)-1)*100, P.eff2=(exp(P.log2)-1)*100, RH.eff2=(exp(RH.log2)-1)*100,
                                   TRH.eff1=(exp(TRH.log1)-1)*100, TRH.eff2=(exp(TRH.log2)-1)*100, RHP.eff1=(exp(RHP.log1)-1)*100,RHP.eff2=(exp(RHP.log2)-1)*100,
                                   Tot.eff1=(exp(Tot.log1)-1)*100, Tot.eff2=(exp(Tot.log2)-1)*100, 
                                   yield0, area), by=.(Crop,Irrig,Level)]

#merge effects of 8 crops to 6 crops and merge irrig/rainfed
sens.sai.pct[Crop=="Tropical Corn", Crop:="Corn"]
sens.sai.pct[Crop=="Tropical Soybean", Crop:="Soybean"]

#per crop effect: merge irrig/rainfed, weight by area and base yield
sens.sai <- sens.sai.pct[, .(T.eff1=weighted.mean(T.eff1, yield0*area, na.rm = T), RD.eff1=weighted.mean(RD.eff1, yield0*area, na.rm = T),RI.eff1=weighted.mean(RI.eff1, yield0*area, na.rm = T),
                             P.eff1=weighted.mean(P.eff1, yield0*area, na.rm = T), RH.eff1=weighted.mean(RH.eff1, yield0*area, na.rm = T),Tot.eff1=weighted.mean(Tot.eff1, yield0*area, na.rm = T), 
                             T.eff2=weighted.mean(T.eff2, yield0*area, na.rm = T), RD.eff2=weighted.mean(RD.eff2, yield0*area, na.rm = T),RI.eff2=weighted.mean(RI.eff2, yield0*area, na.rm = T),
                             P.eff2=weighted.mean(P.eff2, yield0*area, na.rm = T), RH.eff2=weighted.mean(RH.eff2, yield0*area, na.rm = T),
                             TRH.eff1=weighted.mean(TRH.eff1, yield0*area, na.rm = T), TRH.eff2=weighted.mean(TRH.eff2, yield0*area, na.rm = T),
                             #RDRI.eff1=weighted.mean(RDRI.eff1, yield0*area, na.rm = T), RDRI.eff2=weighted.mean(RDRI.eff2, yield0*area, na.rm = T),
                             RHP.eff1=weighted.mean(RHP.eff1, yield0*area, na.rm = T), RHP.eff2=weighted.mean(RHP.eff2, yield0*area, na.rm = T),
                             Tot.eff2=weighted.mean(Tot.eff2, yield0*area, na.rm = T), area=sum(area, na.rm = T)
                             ), by =.(Crop,Level)]

sens.sai #negative insolation effect for all crops per unit changes in RD,RI
#======report these values and the sensitivity figure to answer Reviewer#1 and #3 to show that our model coefficients can reproduce the effects of Pinatubo eruption as in Proctor et al,
#======but the AOD data of Pinatubo is an extreme case, and might need to remove background aerosol effects because +0.15 AOD is much higher than CMIP6 volcanic only AOD for 1992 and 1983


sens.sai1 <- sens.sai[,.(T.eff=T.eff1, RD.eff=RD.eff1, RI.eff=RI.eff1, P.eff=P.eff1, RH.eff=RH.eff1, TRH.eff=TRH.eff1, RHP.eff=RHP.eff1
                         ), by =.(Crop,Level)]
sens.sai2 <- sens.sai[,.(T.eff=T.eff2, RD.eff=RD.eff2, RI.eff=RI.eff2, P.eff=P.eff2, RH.eff=RH.eff2, TRH.eff=TRH.eff2, RHP.eff=RHP.eff2
                         ), by =.(Crop,Level)]
sens.sai1.t <- sens.sai1[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Crop)]
names(sens.sai1.t)[2:5] <- c("Variable", "Mean", "Conf.L",  "Conf.U")
sens.sai2.t <- sens.sai2[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Level", by=.(Crop)]
names(sens.sai2.t)[2:5] <- c("Variable", "Mean", "Conf.L",  "Conf.U")


(crop6 = unique(sens.sai1.t$Crop))
crop6 = c("Corn", "Sugarcane", "Wheat", "Rice", "Soybean", "Cotton")
cropname = c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")
(vars = unique(sens.sai1.t$Variable))
varname1 = c("-1 [°C] T", "-10% RD", "+10% RI", "-10% P", "+5 [%] RH" , "-5 T*RH", "-0.5 RH*P")
varname2 = c("-1sd T", "-1sd RD", "+1sd RI", "-1sd P", "+1sd RH" , "-1sd*sd T*RH", "-1sd*sd RH*P") 
# library(RColorBrewer)
# cbPalette <- brewer.pal(n =12, "Paired") # ,"Dark2")
# varcolor <- c(cbPalette[c(10,7,8,1,2,3,5)], "#000000", "#737373")
# barplot(height=rep(1,length(varcolor)),col=varcolor) #display colors

#====final colors for variables and scenarios
varcolor3 <- c("#E64B35FF","#FF7F00" ,"#FDBF6F", "#A6CEE3", "#1F78B4", "#B2DF8A", "#FB9A99", "#000000", "#737373")
barplot(height=rep(1,length(varcolor3)),col=varcolor3) #display colors


jpeg(paste0(wd1,"/figures/extended/ED_Fig4a.jpg"), units = "in", width = 4.5, height = 3, res = 600)
#svg(paste0(wd1,"/figures/supplementary/Fig.S3a.svg"), width = 4.5, height = 3)

p <- sens.sai1.t  %>% 
  mutate(Crop = factor(Crop, levels=crop6, labels=cropname)) %>%
  mutate(Variable = factor(Variable, levels=vars, labels=varname1)) %>%
  ggplot(aes(x=Crop, y=Mean, fill=Variable)) + 
  #ggplot(aes(x=Variable, y=Mean, fill=Crop))+ # , alpha=Irrig)) +
  geom_bar(stat="identity", width=0.75, position=position_dodge(width=0.75)) + #position=position_dodge(width=400), binwidth=500) +
  geom_errorbar(aes(ymin=Conf.L, ymax=Conf.U), width=.2,
                position=position_dodge(width=.75)) + #smaller width increase the gap between dodge'd bin groups
  #facet_wrap(~Scenario, ncol=1, labeller= label_parsed) +
  scale_x_discrete(cropname) + # in expression (µE m"^-2~d^-1*") #end superscript ^ with *, add space by ~
  scale_y_continuous(expression(Delta~Yield*" (%)"), limits=c(-15, 48), breaks = c(-10, 0, 10, 20, 30, 40)) + 
  scale_fill_manual(name="", values= varcolor3 )+
  guides(fill = guide_legend(label.position = "right", nrow=2, byrow=T, keyheight=0.6, keywidth = 0.6, label.hjust = 0))  #use keyheight/keywidth to adjust legend box size
p + theme_bw()+  
  theme(title = element_text(size = 10)) + #ggtitle(expression('('*bold(a)*') SAI - RCP85')) +
  ggtitle(expression(''*bold(a)*'')) + 
  theme(legend.text = element_text(size=10),  legend.position = "bottom") +  #legend.margin = margin(6, 6, 6, 6),
  theme(axis.title=element_text(size=10), axis.text = element_text(size=10)) +
  theme(axis.title.x =element_blank(), # axis.text.x = element_text(angle=45,hjust=1), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

dev.off()


jpeg(paste0(wd1,"/figures/extended/ED_Fig4b.jpg"), units = "in", width = 4.5, height = 3, res = 600)
#svg(paste0(wd1,"/figures/supplementary/Fig.S3b.svg"), width = 4.5, height = 3)

p <- sens.sai2.t  %>% 
  mutate(Crop = factor(Crop, levels=crop6, labels=cropname)) %>%
  mutate(Variable = factor(Variable, levels=vars, labels=varname2)) %>%
  ggplot(aes(x=Crop, y=Mean, fill=Variable)) + 
  geom_bar(stat="identity", width=0.75, position=position_dodge(width=0.75)) + #position=position_dodge(width=400), binwidth=500) +
  geom_errorbar(aes(ymin=Conf.L, ymax=Conf.U), width=.2,
                position=position_dodge(width=.75)) + #smaller width increase the gap between dodge'd bin groups
  #facet_wrap(~Scenario, ncol=1, labeller= label_parsed) +
  scale_x_discrete(cropname) + # in expression (µE m"^-2~d^-1*") #end superscript ^ with *, add space by ~
  scale_y_continuous(expression(Delta~Yield*" (%)"))+#, limits=c(-15, 26), breaks = c(-10, 0, 10, 20)) + 
  scale_fill_manual(name="", values= varcolor3 )+
  guides(fill = guide_legend(label.position = "right", nrow=2, byrow=T, keyheight=0.6, keywidth = 0.6, label.hjust = 0))  #use keyheight/keywidth to adjust legend box size
p + theme_bw()+  
  theme(title = element_text(size = 10)) + ggtitle(expression(''*bold(b)*'')) + #ggtitle(expression('('*bold(b)*') SAI - RCP85')) + 
  theme(legend.text = element_text(size=10),  legend.position = "bottom") +  #legend.margin = margin(6, 6, 6, 6),
  theme(axis.title=element_text(size=10), axis.text = element_text(size=10)) +
  theme(axis.title.x =element_blank(), # axis.text.x = element_text(angle=45,hjust=1), 
        panel.grid.major = element_blank(), panel.grid.minor = element_blank())

dev.off()







